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1. Introduction 


The atomization and breakup process of liquid fuel jet is of vital importance to 
combustion performance of many practical devices, such as liquid rocket engines, gas-turbine 
combustors and diesel engines. In SSME or STME, for example, non uniform mixing of the 
liquid oxygen/hydrogen propellants in the atomization processes is responsible for the 
injector performance loss and influences the consequent evaporation, gas-gas mixing and 
combustion, and combustion instability. 

Numerical modeling of liquid-jet atomization requires the resolution of the different 
conservation equations that describe the dynamics of the interfaces separating two 
immiscible fluids as well as the droplet/spray dynamics due to the breakup of liquid droplets 
from the liquid-gas interface. Past atomization models [1,2] have relied highly on the simple 
one-dimensional idealized model [3,4] or using empirical formulation based on specific 
injectors. 

During the last decade, general multiphase computational fluid dynamics 
methodologies have matured to the point that several attempts have been made to model the 
detailed atomization process. Reitz [5] has developed an "Blob Injection" atomization model 
based on wave instability analysis. This model has been casted in the KIVA code for diesel 
engine applications^] and also been used by Kim et al. [7,8,9] to study combustion 
instability of a annular liquid combustor model Habiballah et al. [10]. Przekwas and co- 
workers [11] have used direct numerical simulation to study the instability and breakup of 
laminar jets and have device a "jet embedding" technique to couple with flow solvers for 
liquid-engine thrust chamber calculations. Probably the most rigorous numerical modeling of 
atomization process are the works done by Liang and co-workers [12,13,14,15]. 

The ARICC-3D code developed by them combined the volume of fluid (VOF) 
methodology with the spray combustion code KIVA II. The VOF method is used for tracking 
two immiscible fluids and the resulting spray/droplets and combusting gas are solved by 
Lagrangian particle tracking and ALE/ICE (Arbitrary Lagrangian-Eulerian/Implicit 
Continuous Eulerian) finite volume differencing of the reacting flow equations. The ARICC-3D 
enables the simultaneous treatments of three phases: a compressible gaseous mixture, an 
incompressible liquid, and a dispersed droplet spray. Recent applications [15] involve a 
single element injector analysis and a multi-element injector analysis to study the combustion 
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response to a bomb blast with and without baffles. The multi-element computations were 
found very time-consuming. The authors lists several areas where effort has to be made to 
enhance the computational efficiency. One of the improvement which they made recently is to 
recast the VOF formulation into a pressure-based flow solver based on the SIMPLE 
algorithm [16]. 

In the last few years, we have used the state-of-art pressure based method, the 
PISOC algorithm, for calculating chemical reacting flows at all speeds involving spray 
combustion [7,8,9,17,18]. Various physical models including the non-isotropic algebraic 
stress model (ASM), the multi-scale model, equilibrium and finite rate chemistry, turbulent 
modulation due to droplets, group droplet dispersion due to turbulence, droplet coalescence 
and breakup have been incorporated into the current MAST code [19,20]. 

The MAST code uses primitive variables on a non-staggered general curvilinear grid 
system and follows the PISOC algorithm on a time-marching scheme. One predictor- 
muticorrector sequence was formulated within each time step for time-accurate transient 
calculations. The Chakravarthy-Osher's high order scheme [21] was used for the convection 
terms in the governing equations and the conjugate- gradient (CGS) matrix solver was used 
for solving the descretized algebraic equations sequentially for each variables with the 
predictor-corrector sequence. The purpose of this study is to extend this algorithm to involve 
volume-displacement effect encountered in very dense region in the primary atomization 
regions. The fractional volume of fluid (VOF) method will be coupled with the existing 
Eulerian-Lagrangian scheme currently used in the MAST code to resolve three phases: an 
incompressible liquid fuel phase, a compressible gas phase, and a dispersed droplet phase, 
within the calculation domain. To calculate the surface tension effect, we used the continuum 
surface force (CSF) model [22,23]. This model interprets surface tension as a continuous, 
three-dimensional effect across an interface, rather than as a boundary condition on the 
interface. The continuum method eliminates the need for interface reconstruction, and 
simplifies the calculation of surface tension. 

In this study, to verify the tracking of free surfaces between liquid and gas phases and 
to analyze the interfacial phenomena between liquid and gas phases, we assume the gas and 
liquid phases are incompressible. The confined dam broken problem and water sloshing 
problem were carried out. Also, to verify the surface tension force effects, the single droplet 
problem and the jet breakup problem were solved. Detailed formulation and validation will be 
described in the following sections. 


2 



2. Governing Equations and Physical Models 
2-1. Governing Equations 


The mathematical formulations for the three-phases (liquid, gas, and droplet phase) 
flow and combustion processes comprise the Eulerian conservation equations for the liquid 
and gas phases and Lagrangian equations for the fuel droplets. The link between three 
phases is mathematically expressed in terms of the interaction source terms in the governing 
equations [19]. The tracking of the free surface between the liquid and gas phase is 
represented by the VOF method. 

All phases processes are modeled by a system of unsteady, multi-dimensional 
equations. The gas and liquid phases are written in Eulerian coordinates whereas the liquid- 
droplet phase is presented in Lagrangian coordinates. The two-way coupling between the 
two phases is described by the interaction source terms which represent the rates of 
momentum, mas and heat exchange. The detailed equations can be found in ref. 19. The 
current method is intended to predict the motion of the fluid interfaces based on the use of a 
conserved scalar variable transport equation. The conserved scalar is the fractional volume of 
fluid (VOF) cell partitioning function, and the solution of which provides information on the 
position and shape of the interface [24,25]. Through a linear relation, it also determine the 
fluid properties. By defining the fractional volume in a typical control volume cell: 

F= _ (2-D 

Vg + Vl 


where V represents volumes occupied by gas phase (Vg) or liquid phase (Vi) within the 
control volume considered. In the absence of interfacial heat and mass transfer, the function F 
obeys the volume flux conservation equation: 

^ + =L( U jF)= 0 (2-2) 

3t dxj 

For this initial study, by assuming both liquid and gas flows are incompressible, all other 
conservation equations can be expressed in terms of volume of fraction F. These are 
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continuity equation. 


dp d , — . n 

* + S[ ( P Ui) = 0 


(2-3) 


and momentum conservation equation. 




dt dx 




“J 


(2-4) 


where p is the averaged density defined as; 

p= piF + (1 -F) p g 


(2-5) 


and, gi is the body force (gravity), F sv is the volume force for the surface tension effects, and 
k is the turbulent kinetic energy and the viscous stress tensor is 


X ij= ^ 


dll; 


dll; 


dll 


dx; dxj 3 dxj£ 


Mi 


( 2 - 6 ) 


2-2. Physical Models 

To analyze the atomization of liquid rocket engine and spray combustion, it is 
necessary to incorporate the knowledge and concepts of liquid-gas, droplet-gas, and liquid- 
droplet-gas interfacial phenomena and modeling for resolving the liquid volumes displacement 
effects. To resolve the dynamics of the interfaces separating two immisciable fluids, the 
tracking of free surfaces has to be considered. 

In the present study, the ffee-surface tracking procedure to analyze the interfacial 
phenomena between the liquid and gas-phase in atomization and spray combustion is 
represented by the volume of fraction (VOF) method. And also, for the calculation of surface 
tension force, the continuum surface force (CSF) model [22,23] is used. The various physical 
models including the non-isotropic algebra stress model (ASM), the multi-scale model, 
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equilibrium and finite rate chemistry, turbulent modulation due to group droplet dispersion, 
droplet coalescence and breakup can be found in refs. 7,8,9,17 for details. 

a. VOF Method for Tracking the Free Surface 

In fluid dynamics, Lagrangian and Eulerian coordinates have been used to the tracking 
of free surfaces. In the first approach using Lagrangian discrete representation of a fluid with 
free surface, each zone of grid that subdivides the fluid into elements remains identified with 
the same fluid element for all time. In this case, the grid moves with the computed element 
velocities. It has the advantage of circumventing the problem of numerical diffusion across the 
interface. However, the methods become inapplicable whenever the deformation of the 
interface is severe, such as in droplet and liquid jet breakup studies. 

The second method seeks to retain the numerical versatility of a purely Eulerian 
representation. However, the convective flux calculation requires an averaging of the flow 
properties of all fluid elements that fluid themselves in a given mesh cell after some period of 
time. The convective averaging results in a smoothing of all variations in flow quantities, and 
in particular, a smearing of surfaces of discontinuity such as free surfaces. The only way to 
overcome this loss in boundary resolution is to introduce some special treatment that 
recognizes a discontinuity and avoids averaging across it. 

One of the special treatments is the VOF (Volume of Fraction) method. This method 
was developed by the Los Alamos Group [24,25]. This method forms the basis of the SOLA- 
VOF program [25]. The SOLA-VOF solution algorithm has been designated for a wide range 
of applications. It may be applied to problems involving a single fluid having any number of free 
surfaces, or to two immiscible fluids separated by any number of free interfaces. In this 
technique, a function F(x,y,t) is defined whose value is unity at any point occupied by fluid and 
zero elsewhere. In particular, a unity value of F corresponds to a cell full of fluid, whereas a 
zero value indicates that the cell contains no fluid cells with F values between zero and unity 
contain a free surface. 

Let F(x,y,t) stands for a conserved scalar variable of the fluid that define the fractional 
volume in a typical control volume cell (Eq.(2-1)). That is, free surface reconstructed by means 
of a conserved scalar variable F(x,y,t), where 
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F(x,y,t) = 


(2-7) 


1 : in the liquid phase 

0< F(x,y,t) <1 : at the free surface 

0 : in the gas phase 

Consider at a fixed point r© in space, we can obtain F as a function of time and can 
therefore calculate the rate of change dF/dt at this point. From this limited knowledge of F we 
have no way of knowing, how F changes with time if we stay with a particular panicle and 
following it along as it passes through the point at r 0 . To find this Lagrangian rate of change 
we need to relate the value of F at r Q at time to to its value at a neighboring point r 0 + dr at 
time to+dt, where dr= v dt is a small displacement along the flow line passing through the 
point at r 0 . Now F has the value F 0 = F(xo, y 0 , to) at r 0 at time to- When a panicle at this 
point at the time to arrives at the neighboring point at the time to+dt, the function F has the 
value 


F(xo+d^, y 0 +dr|, t+dt) 


( 2 - 8 ) 


dF 


dxl 0 


dF\ 


dylo 


-Fo+|£U*|Stidn + |£u 


at 


(2-9) 


Where, t, and r) represent the displacement in the direction of x and y respectively. The total 
increment in F is therefore 


dF = 





dt 


( 2 - 10 ) 


Hence, the time rate of change of F from the Lagrangian viewpoint is 
dE = | 9§. + |3F\ dti /dF\ 

dt lax/o 3t lay/o dt | at) 


( 2 - 11 ) 


The conserved scalar variable F(x,y,t) is governed by the following transpon equation 


dF dF 

-5- + Ui — 
dt 3xj 


= 0 


( 2 - 12 ) 


where F(x,y,t) at time t=0 has to be given. Although the Eq.(2-12) is obtained by the 
particular point and at a particular time, both position and time and arbitrarily chosen; this 
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Eq.(2-12) constitutes a general kinematical relation that always holds for a fluid. However, in 
this study, by assuming both liquid and gas flows are incompressible. Eq.(2-12) can be written 
in the form of Eq.(2-2) 

When the VOF equation is solved over a computational cell, the changes in F in a cell 
reduce to fluxes of F across the cell faces. As previously noted, to calculate VOF value, we 
need a special numerical technique in computing these flux to reserve the sharp definition of 
free boundaries. Several researchers have previously used variations of this approach for 
tracking material interfaces. In SOLA- VOF [25], donor- acceptor flux approximation was used. 
The basic idea of this method is to use information about F downstream as well as upstream 
of the flux boundary to establish a crude interface shape, and then to use this shape in 
computing the flux. In ref. [26], the Van-Leer method was used to reduce the numerical 
diffusion and to calculate the flux accurately. In this present studies, the higher order 
Chakravarthy-Osher scheme [21] is used to solve the VOF equation. 

b. CSF Model for Surface Tension Force. 

Liquid surfaces are in a state of tension, as though they possessed an elastic skin, 
because fluid molecules at or near the surface experience uneven molecular forces of 
attraction. Since abrupt changes in molecular forces occur when fluid properties change 
discontinuously, surface tension is an inherent characteristic of material interfaces. Surface 
tension results in microscopic, localized "surface force" that exerts itself on fluid elements at 
interfaces in both the normal and tangential directions. Fluid interfacial motion induced by 
surface tension plays a fundamental role in many natural and industrial phenomena [23]. 

The free boundary between the liquid and gas is known as the interfacial region or, 
simply the interface. The interface region is that thin layer surrounding a geometric surface of 
separation, within which the physical properties differ noticeably from those in either of the 
bulk phases. The thickness of this layer is ill-defined because the variation of physical 
properties across it is continuous. To calculate the surface effects at the interfaces between 
two immiscible fluids, previous researcher adopt an approximation in which the interface is 
infinitely thin; that is, they regard the phase boundary as a geometry surface, and assume 
that the properties right up to the interface are unchanged from those of the respective bulk 
phase. An alternative description from an energetic point of view follows from the fact that 
because a liquid molecule at a liquid-gas interface must be attracted to less neighboring 
molecules than one in the interior of the fluid, the attractive energy per molecule at the surface 
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must then be some fraction of that in the interior. The energy of a surface molecule is 
therefore higher than that of one in the bulk liquid, so energy must be expended to move a 
molecule from the interior to the surface. However, since the free energy of the system will 
tend toward a minimum, the surfaces of the liquid phase tend to contact. With a the force per 
unit length tending to contract the surface, we may therefore write that, at constant 
temperature and volume for a given number of moles of system, 



where, G represents the free energy and A represents surface area. The quantity a is called 
the surface tension and is usually given in units of force per unit length. For a liquid-gas 
interface problem, there is an imbalance of intermolecular forces, although smaller, the 
magnitude of the interfacial tension usually lies between the surface tensions of each liquid. 
And, there will be a tendency to curve the interface, as a consequence of which there must be 
a pressure difference across the surface with the highest pressure on the concave side. The 
expression relating this pressure difference to the curvature of the surface is usually referred 
to as the Young-Laplace Equation. From a calculation of the P-V work required to expand the 
curved surface and so change its surface area, it is relatively straightforward to show that 
this equation may be written as 

AP = ° ( Ri' + Rl ) (2 ‘ 14) 

where Rq and R 2 are the radii of curvature of the surface along any two orthogonal tangents ( 
principal radii of curvature), and AP is the difference in fluid pressure across the curved 
surface. Note that the individual contribution of either Ri or R 2 to the pressure difference is 
negative when moving radially outward from the corresponding center of curvature. As Eq.(2- 
14) is written, it is applicable to arbitrarily shaped surfaces where the radii of curvature may 
change spatially. 

However, Eq.(2-14) has suffered from difficulties in modeling topological complex 
interface having surface tension. In this smdy, surface tension at free surface is modeled with 
a localized volume force prescribed by the recent CSF (continuum surface force) model. In 
CSF model, instead of a surface tensile force or a surface pressure boundary condition applied 
at a discontinuity, a volume force due to surface tension acts on fluid elements lying within 


8 



finite thickness transition regions replacing the discontinuities. CSF formulation makes use of 
the fact that numerical models of discontinuities in finite volume and finite difference scheme 
are really continuous transitions within which the fluid properties vary smoothly from one fluid 
to another. The volume force in CSF model is easily calculated by taking first and second 
order spatial derivatives of the characteristic data. In the case of theVOF method, it is the 
VOF function F. At each point within the free surface transition region, a cell-centered value 
F sv is defined which is proportional to the curvature K of the constant VOF surface at that 
point. 


Surface tension modeled with the continuum method eliminates the need for interface 
reconstruction, so restriction on the number, complexity, or dynamic evolution of interfaces 
having surface tension are not imposed. 

Surface Tension Force 

The surface stress boundary condition at an interface between two fluid is 

(Pi - P 2 +G K)llj =(Tiin - X 2 ik>lk (2-15) 

dXj 

where a is the fluid surface tension coefficient. Pi and P2 are the pressure in fluid(l and 2). T 
is the viscous stress term. K is the free surface mean curvature. In this study, for the accurate 
modeling of the normal boundary condition for interface, we assumed that the viscosity at the 
interface is neglected and the surface tension coefficient is constant. Therefore, the fluid 
pressure jump across an interface under surface tension is 

P S = P 2 - Pi = a <x) (2-16) 

where P$ is the surface pressure. Surface pressure is therefore proportional to the curvature 
of the interface. Surface tension contributes a surface pressure that is the normal force per 
unit interfacial area. Therefore, the surface force per unit interfacial area can be written as 

Fsa(x) = a K(x)n (2-17) 

where n is the normal vector at the free surface. In the CSF model, the surface tension is 
reformulated as a volume force Fsv satisfying 
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lim 

h-»« 


AS 


(2-18) 


F sv (x)d 3 x = 


Fsa(x)ds 


and h is a length comparable to the resolution afforded by a computational mesh with spacing 
dx.(Fig.l) The area integral is over the portion AS of the surface lying within the small 
volume of integration AV . The finite difference approximation in MAST-VOF replace free 

surface discontinuities with finite thickness transition regions within which the fluid 
properties vary smoothly from fluid to gas over a distance of 0(h). 


F sv (x)= OK(x)— (2-19) 

[F] 

where F is the fluid characteristics, equal with the VOF value in MAST-VOF. When F= F, 
the volume force is computed accurately for any two fluids meeting at the interface. In 
particular, the two fluids could have equal densities. The mean free surface curvature K 
given by [22,23] 

k(x) = -Vn = i[(S- V) fn| -(Vn)] (2-20) 


where the unit vector 



( 2 - 21 ) 


is derived from a normal vector n 


n = VF (2-22) 

that is the gradient of VOF data. The volume force, F S v(x) has the following properties; 

(1). The volume force in the transition region in Fig.l, where there characteristic 
varies smoothly form fluid 1 to fluid 2, is designated to simulate the surface 
pressure on the interface between the fluids. Thus, the line integral of F S v(x) 
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across the transition region is approximately equal to the conventional surface 
pressure: 



(2). In the limit that the width of the transition region in a direction normal to 
the interface goes to zero, the volume force becomes the conventional surface 
pressure. 


Modeling surface tension requires some special consideration, since the effects of 
surface tension should be confined to the neighborhood of the interface. To simplify the 
application of boundary conditions and to localize the domain of dependence of the volume 
force, an approximation with compact support is sought. To maintain the integrity of the 
transition region, the volume force should not change sign along the radius of curvature. 

Because the contribution to the surface tension force come from the small portion of 
the computation mesh in the neighborhood of the interface, difficulty in formulating sufficiently 
accurate finite difference expression might be expected. It turns out that low-order 
approximations may be used, provided one begins with a form of the volume force that 
emphasizes the region of maximum gradient. This allows one to apply boundary conditions 
with no more difficulty than with other terms, such as pressure. In the CSF model for surface 
tension, a surface force is formulated to model numerically surface tension effects at fluid 
interface having finite thickness. The method is basically suited for Eulerian interfaces that 
are in general aligned with the computational grid. It can alleviate previous topological 
constraints on modeling interfaces having surface tension without sacrificing accuracy. 

The CSF model has been validated on the single droplet problem at the equilibrium 
state. This problem will be described in the next section. 
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Computational grid 
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Fig.l Schematic Diagram for Interface region 


3. Finite Difference Equation 

3-1. Finite Difference Equation for Transport Equations 


The control volume method was used to derive the difference equations with all 
variables located at each cross point of a grid mesh. The grid system is shown in the Fig. 2, 
where the dash lines are the boundaries of control volume with east, west, north and south 
faces and P, E, W, N, S are grid main points. The transport equations are discretized by the 
Euler implicit difference scheme. The governing equations can be expressed in difference form 
for each grid point as 

Continuity equation: 

-J- [ p«' - p" } + Ai (pur 1 = 0 (3-1) 

At 

Momentum equation: 

j- [ (p Ur 1 - (p U) n ] = H(U” +1 ) - Ai P"“ + S u + S„ + S M (3-2) 
At 

VOF equation 

_1_ [ p +1 - p ] + A* (Ui Ff +1 = 0 (3-3) 

At 

In the above equations, the operators H denotes the finite difference representation of the 
spatial convective and diffusive fluxes of velocity Ui, F stands for VOF value. The operator 
A i represents the first order Euler finite difference of 3/dx,. The source terms, S u , contains all 
other terms except the convective and diffusive term for each variables. S st stands for the 
surface force due to surface tension. Sbf represents body force due to gravity. 

In order for the solution procedure of the finite difference equations to be stable, the 
simplest way without losing accuracy is to separate the diagonal elements of the operators H’ 
and to shift them to the left-hand side of the equations. The concentration is focused on the 
momentum equation: 
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where. 


H (UO = H (UO - Ap Ui 


(3-4) 


H operator of convection and diffusion term at neighbor points 
surrounding the main point p 

Ap coefficient of the diagonal element of the operator H. 

3-2. Finite Difference Equation for VOF Equation 

In Eulerian representation, the convective flux calculation usually requires an 
averaging of the flow properties of all fluid elements in given mesh cells. This convective 
averaging results in a smoothing of all variations in flow quantities, and in particular, a 
smearing of surfaces of discontinuity such as free surfaces. To calculate the interface fluxes 
for the advection term, the upwind scheme is too dissipative, while the central and Lax- 
Wendroff schemes are too dispersive in the vicinity of discontinuities. To overcome this loss 
of accuracy in boundary resolution and to handle sharp interfaces between liquid and gas 
phases, a high accuracy scheme is implemented. This high accuracy scheme must have the 
ability to generate numerical algorithm which allow a high resolution of discontinuities, such 
as shock waves and contact discontinuities, without oscillation. In this study, to prevent the 
generation of numerical oscillations and to calculate the interface flues exactly at the 
interfaces between the liquid and gas, the Chakravarthy-Osher scheme [21] is used. This 
scheme is a kind of TVD (Total Variation Diminishing) scheme, whereby the variation of the 
numerical solution is controlled in a non-linear way, such that it forbids the appearance of any 
new extremum.. This scheme can be defined essentially in terms of one parameter. By 
various choices of this parameter, one can obtain schemes with a wide range of accuracy 
including high accuracy (low truncation error) second-order schemes, the conventional 
second-order accurate upwind TVD scheme and even a third-order accurate TVD scheme. 
This scheme can easily apply to scalar equations, to systems of equations, to arbitrary 
curvilinear coordinate systems, and to general control volumes. 

The transport equation for the conserved scalar variable F can be transformed to a 
general curvilinear coordinate system. From Eq. (2-2), this equation can be written; 

J S + | <UF| + li (VF, = 0 < 3 - 5) 
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Where, U and V are the contra-variant velocity in the transformed coordinates. 
Integrate with respect to the control volume in the grid system Fig.(l). 


J P *' A ; F* + |UF|e - |UF] W + |VF] n - |VFt = 0 


In the above equation, the interface fluxes are calculated by using the Chakravarthy-Osher 
scheme. The F flux term in the direction of \ is obtained as follow; 



_fi+l/2j-fi-l/2j 


(3-7) 


then the control-volume interface flux is computed according to; 


r i+l/2j 


= h 


i+l/2j 


d+O) 




(1 + 0 ) 


4 ^ 

, (l-®) Hf f d-®) 

+ — 1 — dt i-Wj - — 


dfr-wj 


df M/2j 


In the above, h represents a first-order numerical flux, and can be expressed 
hi+l/2j = Uj+l/2j Fij + Uf-l/2j Fi+1 j 

where 

l£i/2j = 0.5 ( U i+1/2J ± |Ui„/2j ) 

The flux-limited values of df are computed as follows: 

dff + i/2,j = Uf + i/2j ( Ff +1J - FPj ) 

= Ut+i/2j ( F?j - F?-ij ) 
dfi + i/2j = Ul+ 1/2 j ( F? + ij - F^j) 
dfi+3/2J = U*i+l/2j ( F?+2j - F?+l j ) 


(3-8) 


(3-9) 


(3-10) 


(3-11) 

(3-12) 

(3-13) 

(3-14) 
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Also, the west, north, and south face fluxes can be obtained by a similar way. 

The spatial accuracy of the scheme is controlled by the parameter 0 , which may take the 
following values: 


0 

= -1 

Fully upwind scheme 

0 

= -1/3 

No name scheme 

0 

= 0 

Fromm scheme 

0 

= 1/3 

Third order upwind scheme 

0 

= 1/2 

Low truncation error second order scheme 

0 

= 1 

Central difference scheme 


3-3. Finite Difference Equation for Surface Tension Force 


The CSF model for surface tension is implemented by placing the normal at grid main 
point and the curvature K at cell centers as shown Fig.l. The shading part is the control 
volume for the computational domain. 

From Eq.(2-19) and Eq.(2-20), the vector and the mean curvature also can be 
expressed 


- 0F dF 
n "8x ei ej 


K(X) = L 

M 




( n * n y\ 

|3n x 

♦ft 

f [dayf 

iff J 

1 dy ax J 

\PI. 

1 Uy) 


(X ILxdF) + 

\r a 3y\ dyl dx\dxl 


(3-15) 


(3-16) 


And the derivatives in above equations are transformed to a general form based on a non- 
orthogonal coordinate system (^ , Tj) . The first derivatives are as follows; 



(3-18) 


8F _ v 9F 3F 

aT^ + T1 y^ 

The second derivatives are as follows; 



The cross derivatives are as follows; 



And the second derivative for the cylindrical coordinate is as follow; 



(3-19) 

(3-20) 

(3-21) 

(3-22) 

(3-23) 



and, the superscript OC on the radius r is a constant equal to 1 in the cylindrical coordinate 
and 0 in the Cartesian coordinate. The finite difference equations for above derivatives are 
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based on the control volume in the grid system (Fig.2). The finite difference equations are as 
follows; 


dF 

a? 

i 

+ 

II 

<u 

Fij 

(3-26) 

dF 

= Fy ' 

w 

Fi-ij 

(3-27) 

dF 

dn 

= F ij+l - 
n 

Fy 

(3-28) 

dF 

dn 

= Fy - 

s 

Fij-i 

(3-29) 


therefore, the first derivatives can be expressed 


[ n x]p ~ 


dF 


dx 


JP 


~2 ( ^x,e 


dF 




+ ^x, 




dF 


T lx4i 


L^Jw L^lJ 


dF 


In 


+ 7 lx,s 


[ n y]p ~ 


® = i (i 

ayjp 2 




dF 


L^Jw y m 


^ydi 


dF 


+ 


Jn 




dF 

^lj 


dF 


(3-30) 


(3-31) 


and, the second order derivatives are; 


dn x 


= [n Je - [n x ]p 


* 

dn x 


= [n x ]p -[n x ]w 


Iw 


^lh 


— [Hx]n ' [^x]p 


(3-32) 


(3-33) 


(3-34) 
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(3-36) 


Similarly, the second order derivative in the direction of y can be obtained; 



(3-39) 

These discretized equations are then substituted into Eq. (3-2) for calculating body 
forces based on the CSF model. 
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4. Numerical procedure 

In this study, the basic algorithm for flow solver is the PISOC algorithm [20], which is 
based on the operator-splitting technique originally proposed by Issa [28], that is, the 
pressure implicit with splitting of operators (PISO) algorithm. PISO algorithm is a non- 
iterative method for handling the coupling of the implicitly discretized time-dependent single 
phase fluid flow and heat transfer. The PISO algorithm has been extended to include a 
stochastic Lagrangian panicle tracking scheme for computing time-dependent gas-droplet 
flows [28], Since for this study both gas and liquid phase are assumed incompressible, the 
variation of density is solely accounted for by the VOF information. By substituting the 
volume flux conservation equation (2-2) into the continuity equation (2-3). It can be shown 
the continuity equation retains the original incompressible form. Thus the same procedure 
described in [28] can be used here with an addition of a F-equation corrector step at the end 
of the last momentum corrector step. Special attention was paid to the solution of the F- 
equation. 

4-1. PISO- VOF Algorithm 

PISO algorithm is a non-iterative method for handling the coupling of the implicitly 
time dependent fluid flow equations. This algorithm is based on the use of pressure and 
velocity as dependent variables and is hence applicable to both the compressible and 
incompressible versions of the transport equations. The main feature of the technique of the 
splitting of the solution process into a series of steps whereby operations on pressure are 
decoupled from those velocity at each step, with the split sets of equations being amenable to 
solution by standard techniques. At each time-step, the procedure yields solutions which 
approximate the exact solution of the difference equations. The accuracy of this splitting 
procedure is assessed for a linearised form of the discretized equations, and the analysis 
indicates that the solution yielded by it differs from the exact solution of the difference 
equations by terms proportional to the powers of the time-step size. Detailed PISO algorithm 
and the PISOC algorithm can be found in Ref. [20,28] for dispersed multiphase reacting flows. 
The implementation of the VOF is summarized in the following; 


Predictor Step 



The equation of momentum is solved in this step implicitly, using old time pressure 
and density, as 

( At +Ap) ^ = H(U ' ) ' p Aipn + lf +Su i (4-1) 

The solution in this step yields U*. velocity field. But, in this step, U* velocity field will not in 
general satisfy the zero-divergence condition. 


The equation of VOF is solved in this step also implicitly, using and old time VOF 
value F n , as 

( At + Ap) F * = G(F *) + XT < 4 -2> 

where, G term is convective term for VOF. 


First Corrector Step 


In this step, the momentum equation can be written in the explicit corrector form 

( ^+ a p ) u*-= H ( u r)-iA lP *+^ + s Ui ( 4-3) 

which, by subtracting predictor equation from it, and then combining the continuity equation 
and equation of state, the pressure increment equation can be obtained; 

ur=u* - ia, fp-.p") ( 4 .4) 

[ A i (lD u A i )](p , -p'>) = AiUf (4 . 5) 

This equation yields P* field. Form Eq.(4-3), U- * can be obtained. The equation of VOF is 
solved in this step also explicitly, using U* * , F*. and old time VOF value F n , as follows; 

(^ + A p )F" = G<F-) + ^ ( 4.6) 

Second Corrector Step 
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For this step, the momentum corrector equation is 


( -L+ Ap ) UT = H( U**) - -U p"+ -^ + S u . (4-7) 

At P At 

Also, in this step, the increment equation for momentum can be obtained as like the First 
corrector step; 


(-L+ApHUT-UT) 

At 

= H(UT)-H(U*)-lMp"-p*) (4-8) 

Denote 

Du = ( aT +Ap)1 (4 ~ 9) 

Thus, 

U?~ = U" +D U [ ( H( Uf) - H( U* ) ) - r Aj ( p” - p* ) ] (4-10) 


Substitute U *** into the continuity equation, we obtain 

[ Ai ( -t D„ Ai ) ] ( p” - P* ) 

= D„ Ai[ H ( U**) - H ( U” ) ] (4-11) 

The above equation yields P**, and then the other variable, U , p** can be obtained. And, 
also the VOF value is obtained as follows; 

( _L + Ap) F*** = G(F**) + -F- (4-12) 

At At 
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4-2. The Construction of Free Surface 


In this study, it is assumed that the boundary can be constructed by a straight line 
cutting through the cell By determining the slope of this line, it can then be moved across the 
cell to a position that intersects the known amount of F volume in the cell and neighborhood 
cells. 


To determine the free surface position in the cell, first calculate the gradient of volume 
of fraction in the cell which has the value of volume of fraction greater than 0.5, that is, VF 


VF = — 

ax 


®i + 


1J 


3F 

dyy 


(4-13) 


and then, calculate the resultant of VF ; 



3F 2 
+ (|^) 
dy i,j 


(4-14) 


The x-position and y-position of the free surface can approximately obtained from Eq.(4-13) 
and Eq.(4-14) ; 


X sf = Xy + ( 


0.5 - Fij ) ( 


dF/dx y 

~VF| iJ 


(4-15) 


Y S f= yy + ( 0.5- Fg) ( 


aF/9y 

VF| Xj 


(4-16) 


By the similar way, the x and y positions of the neighborhood cells (i-l,j),(i,j-l) and (i+l,j) 
are obtained. And then , if Fjj is greater than 0.5 , Fi+i j is also greater than 0.5, and Fjj+i is 
less than 0.5, connect the point (X s f , Y s f ) of the left-side cell (i-l,j) and the point (X s f , Y s f ) 
of the left-side cell (i,j). In this case, The free surface is constructed horizontally, if Fj j is 
greater than 0.5 , Fi+ij is also less than 0.5, and Fy+i is greater than 0.5, connect the point 
(X s f , Y s f ) of the top-side cell (i,j) and the point (X s f , Y s f ) of the bottom-side cell (i,j-l). In 
this case. The free surface is constructed vertically. 
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5. Validation 


To validate the MAST-VOF, we solved dam broken problem, water sloshing problem, 
single droplet, and liquid column breakup problem. All problems were solved simultaneously 
liquid and gas phase. 

5-1. Dam Broken Problem 


In this calculation, a rectangular column of water in hydrostatic equilibrium is confined 
the wall as shown Fig. 3. At the beginning of the process the right wall is suddenly removed, 
the water column therefore starts to collapse, under the influence of gravity and to form an 
advancing water wave. 



Fig.3 Schematic Diagram for Dam Broken Problem 

This is a good test problem for tracking the free surface because it has simple 
boundary conditions and a simple initial configuration. The appearance of both a vertical and 
horizontal free surface, however, provides a check on the capability of MAST-VOF to treat 
free surface that are not single valued with respect to X-Y coordinate. The physical 
properties are listed in Table 1. 
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Table 1. Physical properties for dam broken problem 


o Liquid physical properties 

density : 1000.0 Kg/m^ 

viscosity : 1.0 E-3 Kg m/sec 
o Gas physical properties 

density : 1.0 Kg/m^ 

viscosity : 1.8 E- 5 Kg m/sec 
o Gravity : 0.01 m/sec^ 

At the start of the calculation, VOF value F has the unit in the region occupied by the 
water, and zero in the gas phase. Also, the velocity components u and v are zero at the flow 
domain boundary. For the calculation, a grid consisting of 41x41 uniformly spaced cells has 
been chosen. The numerical results are shown in Fig. 4. At the initial state, there is no gas 
flow. When the dam collapse, the gas velocity was induced due to the momentum transfer 
from the liquid phase (water). Gas started to move from left to right wall. As the water front 
moves from left to right, the gas velocity increased and the recirculating flow also induced. 
Numerical results from the present method in terms of VOF values are in good agreement 
with those in ref. [26]. Also, the computed results of water column height H and wave front 
length Z are compared with the experimental data of Martin & Moyce [29] as shown in Table 
2. From this comparison, it is seen that the present mathematical model is capable of 
tracking the free surface between the liquid and gas phase. 

Table 2 Comparison of numerical results to experimental data 



MAST-VOF 

Experimental Data 

Time(sec) s N. 

Z(cm) 

H(cm) 

Z(cm) 

H(cm) 

0.0 

1.0 

2.0 

1.0 

2.0 

0.9 

1.514 

1.68 

1.5 

1.7 

1.4 

2.17 

1.37 

2.17 

1.37 

2.0 

3.28 

1.062 

3.1 

1.1 
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Fig. 4 Computational results for dam broken problem 
at time = 0.0 sec 
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5-2 . Water Sloshing Problem 

The second test calculation is that of liquid water sloshing in two-dimensional 
rectangular tank as shown Fig. 5. The tank is left-half filled with liquid water and right-half 
filled with gas. At the initial state, the pressure is uniform at the value Po and the velocity is 
zero. And also, VOF value at the left-half side is equal to unity and the right-half side is zero. 
And the physical properties is listed in Table 3. 

The numerical results are shown in Fig. 6. In this case, no analytic solution is 
available. However, the numerical results clearly show the qualitative sloshing behavior 



Fig. 5 Schematic diagram for water sloshing problem 


which would be expected intuitively. At the beginning state, the water front starts to collapse 
due to the gravity force. As the water front starts to move from left to right, the gas velocity 
was induced due to the momentum transfer from the liquid phase (water). And also the 
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o Liquid physical properties 

density : 960.0 Kg/m^ 

viscosity : 1.0 E-3 Kg m/sec 
o Gas physical properties 

density : 0.585 Kg/m^ 

viscosity : 1.8 E-5 Kg m/sec 
o Gravity : 9.8 m/sec^ 

o Initial Pressure : 1.0 bar 

recirculating flow was induced. The numerical results are compared with ref. [30] with close 
agreement qualitatively. 
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Water Sloshing Problem 
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Ttaie= 0.600000 tyde= 1200000 


Fig. 6 Velocity vector and free surface configuration computed 
for water sloshing problem 
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Fig. 6 Continue 
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Fig. 6 Continue 
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Fig. 6 Continue 
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Fig. 6 Continue 
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To test the surface tension force, the equilibrium droplet problem was chosen. In this 
problem, the viscous, gravitational, or other external forces were not considered. Only the 
surface tension effect is considered. Therefore, this surface tension causes a static liquid drop 
to become spherical. 

For numerical solution, the Cartesian geometry using a two-dimensional 0.06 x 0.06 
(m) computational domain is chosen as shown in Fig. 7. Also, a regular, orthogonal grid with 
uniform mesh is used The physical properties are listed in Table 4 and the radius of droplet is 
0.02 m. This droplet is located at the center on computational domain (0.03m x 0.03m). 



Fig. 7 Schematic diagram for equilibrium droplet problem 
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Table 4. The physical properties for the equilibrium droplet problem 


o Liquid physical properties 

density : 1000.0 

surface tension coefficient 
initial radius : 0.02 

o Gas physical properties 

density : 1.0 

o Gravity : 0.0 


Kg/m 3 

: 0.02361 N/m 
m 

Kg/m 3 
m/sec ^ 


Fig. 8 shows the numerical results, the gradient of F has its greatest magnitude in 
the transition region and falls to zero outside the transition region. It means that the value of 
mean curvature are local to the interfaces between two fluids. To calculate the surface 
tension, the effects of surface tension should be confined to the neighborhood of the interface. 
Because the surface tension is the boundary phenomena. In Fig. 8, it is seen that the velocity 
vectors inside the droplet are in equilibrium. Therefore, the static liquid droplet keeps 
spherical. Also, the surface force vector and the contour of the density are shown in Fig. 8. 
The results were obtained without any smoothing treatment as described in ref. [23]. 
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5-4. Liquid Column Instability Problem 

This problem concerns the instability of a cylindrical column of fluid and a perturbation 
by radial displacements. A cylindrical liquid column 1.4 cm long by 0.1 cm radius was covered 
by a 1.4 cm long and 0.3 cm calculation domain with a 232 x 101 grid. Liquid to background 
fluid density ratio is approximately 2 and surface tension coefficient 59.0 dynes/cm is 
modeled. 

The interfaces between two immiscible is initially perturbed by a cosine function with 
an amplitude of 0.01 Ro (liquid column radius). For large amplitude displacements, the non- 
linearity takes effect and the liquid column deforms and pinch-off occurs to form small satellite 
drops between the large drops. Details of the satellite drop formations are very sensitive to 
the grid system and topological reconstructions of VOF data into the breaking up of liquid 
bulges. A series of computational results of the free surface and velocity filed are shown Fig. 
9 depicting the deformation of the liquid column into droplets. This computational results are 
compared to experimental results of [31] shown in Fig. 10. 

Next, a cylindrical liquid column 1 cm long by 0.05 cm radius was covered by a 1 cm 
long and 0.15 cm calculation domain with a 81 x 31 grid. Liquid to gas density ratio is 
approximately 1000 and surface tension coefficient 72.8 dynes/cm is modeled. The gas-liquid 
interface is initially perturbed by a cosine function with an amplitude of 0.001 a (liquid column 
radius). The normalized perturbation growth plotted in Fig. 11 vs. non-dimensional time is 
compared with the linear theory of Rayleigh. The current results compared favorably with the 
linear theory as well as with other results by SOLA-VOF [25] and ARICC-ST code [13]. 
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Fig. 9 Velocity vector and free surface configuration computed 
for liquid jet instability problem 
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Fig. 9 Continue 
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time = L300M0 cycle = 260 time = L300*10 
















9 


Fig. 10 Experimental results of water jet breakup with a = 59 dynes/cm. 
Time interval between each photograph is 1.25 x 10 3 sec. 

(time increases from left to right, top to bottom) 
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Fig. 1 1 Comparison of numerical results to linear theory for the liquid column 
instability problem ( : linear theory, • : numerical results) 




6. Conclusion and Recommendation 


The development of a combined Eulerian-VOF-Lagrangian method has been initiated 
based on a pressure-based transient gas-droplet solution procedure [28]. This study 
concentrated on implementing the VOF method and CSF model into the solver for sharp 
interface tracking and efficient handling of the surface tension force. Validation studies 
indicated that the current methodology can readily be extended to incorporate other physical 
submodel such as compressibility effects in the gas phase, spray tracking and turbulence for 
atomization process simulation. 

The future work for the next phase studies should include: 

1. Determine the onset of primary breakup using the VOF data and develop an 
criterion coupling with Lagrangian droplet tracking. 

2. Interfacial transport phenomena including heat effects. 

3. Interfacial transport due to turbulence effects. 

4. Extension to three-dimensional multi-jet simulations. This will involve 
extending the VOF method to couple with multi-zone unstructure grid 
methods. 
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